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Abstract 

The finite element methods (FEM) are important techniques in engi¬ 
neering for solving partial differential equations, but they depend heavily 
on element shape quality for stability and good performance. In this pa¬ 
per, we introduce the Adaptive Extended Stencil Finite Element Method 
(AES-FEM) as a means for overcoming this dependence on element shape 
quality. Our method replaces the traditional basis functions with a set 
of generalized Lagrange polynomial (GLP) basis functions, which we con¬ 
struct using local weighted least-squares approximations. The method 
preserves the theoretical framework of FEM, and allows imposing essen¬ 
tial boundary conditions and integrating the stiffness matrix in the same 
way as the classical FEM. In addition, AES-FEM can use higher-degree 
polynomial basis functions than the classical FEM, while virtually pre¬ 
serving the sparsity pattern of the stiffness matrix. We describe the for¬ 
mulation and implementation of AES-FEM, and analyze its consistency 
and stability. We present numerical experiments in both 2D and 3D for 
the Poisson equation and a time-independent convection-diffusion equa¬ 
tion. The numerical results demonstrate that AES-FEM is more accurate 
than linear FEM, is also more efficient than linear FEM in terms of error 
versus runtime, and enables much better stability and faster convergence 
of iterative solvers than linear FEM over poor-quality meshes. 

Key Words: finite element methods; mesh quality; weighted least squares; 
partition of unity; accuracy; stability 


1 Introduction 

The finite element methods (FEM) are arguably one of the most important 
numerical tools for solving partial differential equations (PDE) over complex 


1 


domains in engineering. They account for an overwhelming majority of the 
commercial and research code for modeling and simulations, and there is a vast 
amount of theoretical work to provide a rigorous foundation; see e.g., [T]. 

Despite their apparent success in many applications, classical hnite element 
methods have a very fundamental limitation: they are dependent on element 
shape quality. This is especially true for elliptic and parabolic problems, for 
which the resulting linear system is often ill-conditioned if a mesh contains a 
few “bad” elements. This can lead to very slow convergence of iterative solvers 
and sometimes even a loss of accuracy. Because of this, researchers and users 
of FEM often spend a tremendous amount of time and computing power to 
generate and maintain meshes, trying to fix that one last bad element. This has 
spurred much successful research in meshing, such as Delaunay triangulation 
11 [3], advancing front [4], octree-based methods [5], etc. However, the meshing 
problem continues to become more and more challenging as applications become 
more and more sophisticated and demanding, and also with the increased use 
of higher-order methods, such as spectral element methods [5], discontinuous 
Galerkin methods HIE] , etc. 

The FEM community has long considered this dependency on element qual¬ 
ity as a critical issue, and the community has been actively searching for al¬ 
ternative methods to mitigate the issue for decades. Examples of such alter¬ 
native methods include the diffuse element or element-free Galerkin methods 
unni, least-squares EEM m, generalized or meshless finite different methods 
miiiiiiiisiiis], generalized or extended EEM HZllIHl, and partition-of-unity 
FEM [19] . To reduce the dependency on mesh quality, these methods avoid the 
use of the piecewise-polynomial Lagrange basis functions found in the classical 
FEM. However, they also lose some advantages of the classical EEM. In partic¬ 
ular, for the generalized finite different methods, the strong form instead of the 
weak form of the PDE must be used. The partition-of-unity FEM and other 
similar generalizations often incur complexities in terms of imposing essential 
boundary conditions and/or integrating the stiffness matrix nnj. Therefore, 
it remains an open problem to develop a numerical method that overcomes 
the element-quality dependence, while preserving the theoretical framework of 
FEM, without complicating the imposition of boundary conditions and numer¬ 
ical integration. 

This paper introduces a new method, called the Adaptive Extended Sten¬ 
cil Finite Element Method (AES-FEM) (pronounced as ace-F-E-M), to address 
this open problem. Similar to some of the aforementioned methods, the AES- 
FEM replaces the piecewise-polynomial Lagrange basis functions in the classical 
FEM with alternative basis functions. Different from those methods, our basis 
functions are partition-of-unity polynomial basis functions, constructed based 
on local weighted least squares approximations, over an adaptively selected sten¬ 
cil to ensure stability. We refer to these basis functions as generalized Lagrange 
polynomial (GLP) basis functions. Another difference of AES-FEM from most 
other generalizations of FEM is that AES-FEM preserves the traditional hnite 
element shape functions as the weight functions (a.k.a. test functions) in the 
weak form, to preserve the compact support of integration and the weak form 
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after integration by parts. This combination of the basis and weight functions 
enables AES-FEM to overcome the element-quality dependence, while preserv¬ 
ing the theoretical framework of FEM, without any complication in imposing es¬ 
sential boundary conditions or integrating the stiffness matrix. In addition, the 
resulting stiffness matrix of AES-FEM has virtually the same sparsity pattern 
as that of the classical FEM, while allowing the use of higher-degree polynomials 
and hence significantly improved accuracy. 

As a general method, AES-FEM allows polynomial basis functions of ar¬ 
bitrary degrees. In this paper, we focus on the use of quadratic polynomial 
basis functions, for which the stiffness matrix has virtually the same sparsity 
pattern as the classical FEM with linear basis functions. However, AES-FEM is 
based on the more general weighted-residual formulation instead of the Galerkin 
formulation, and hence the resulting system is nonsymmetric, which is more ex¬ 
pensive to solve than for symmetric matrices. In addition, it is more expensive 
to construct the basis functions of AES-FEM than to use the standard basis 
functions in FEM. Therefore, AES-FEM is conceivably less efficient than FEM 
for a given mesh. However, as we will demonstrate in our experimental results, 
AES-FEM is significantly more accurate than FEM on a given mesh due to 
its use of higher-degree basis functions, and it is also more efficient than FEM 
in terms of error versus runtime. Most impotently, AES-FEM enables better 
stability and faster convergence of iterative solvers than FEM over poor-quality 
meshes. 

The reminder of the paper is organized as follows. In Section [21 we present 
some background information and recent developments of related methods. In 
Section |3l we formulate FEM from a weighted residual perspective, describe 
the construction of generalized Lagrange polynomial basis functions based on 
weighted least squares approximations, and then introduce AES-FEM. In Sec¬ 
tion |31 we analyze the consistency and stability of AES-FEM. In Section [SI we 
discuss some implementation details including the utilized data structure and 
the applicable algorithms. In SectionlHl we present the results of some numerical 
experiments with our approach. Finally, Section [7] concludes the paper with a 
discussion. 


2 Background and Related Work 

In this section, we review some background information and some methods 
closely related to our proposed approach, including the diffuse element method 
and element free Galerkin method, some modern variants or generalizations of 
FEM, as well as the generalized finite different methods. 

2.1 Diffuse Element Method and Element Free Galerkin 
Method 

Various alternatives of finite element methods have been proposed in the lit¬ 
erature to mitigate the mesh-quality dependence. One of the examples is the 
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diffuse element method (DEM) [TD], proposed by Nayroles, Touzot, and Villon 
in 1992. Similar to AES-FEM, DEM constructs local approximations to an un¬ 
known function based on a local weighted least squares fitting at each node. 
However, unlike AES-FEM, the DEM is based on the Galerkin formulation, 
which requires the shape functions to have a compact support for efficiency. To 
this end, DEM relies on a weight function that vanishes at a certain distance 
from a node, in a manner similar to the moving least squares fittings |20) . The 
accuracy and efficiency of numerical integration in DEM depends on the par¬ 
ticular choice of the weight function. In contrast, based on a weighted-residual 
formulation, AES-FEM enforces the compact support of the weak form with 
the weight functions, so the shape function does not need to have a compact 
support. 

Another approach, which is closely related to DEM, is the element-free 
Galerkin method (EFGM) [3], proposed by Belytschko, Lu, and Gu in 1994. 
As a Galerkin method, EFGM also requires a compact support of its shape 
functions, which serve as both the trial functions and weight functions. Similar 
to DEM, EFGM constructs the shape functions based on moving least squares, 
for which the weight function plays a critical role in terms of accuracy and ef- 
hciency. However, depending on the weight functions, the shape functions in 
EFGM may not be polynomials. It requires special quadrature rules with more 
quadrature points than those of standard FEM [9], and it also requires evaluat¬ 
ing a moving least squares fitting at each quadrature point. In addition, EFGM 
requires the use of Lagrange multipliers for essential boundary conditions. In 
contrast, AES-FEM can utilize the same treatments of boundary conditions and 
numerical integration as the standard FEM. 

2.2 Other Variants and Generalizations of FEM 

Besides DEM and EFGM, various other generalizations of FEM have been de¬ 
veloped in recent years. Some of the most notable ones include the generalized 
finite element method (GFEM) and extended finite element method (XFEM) 
[T7l[2Tl[2l[23l[T8]. These methods also alleviate the dependence on the mesh, by 
introducing enrichment functions to replace the standard FEM basis functions 
in regions with discontinuous solutions, such as along cracks. These enrichment 
functions in general are not polynomials, so special quadrature rules may be 
needed for efficiency. Away from the discontinuities, GFEM and XFEM rely 
on the standard FEM discretizations, so good mesh quality is still required in 
general. 

The GFEM and XFEM may be viewed as special cases of the partition of 
unity method (PUM) [19], which is a general framework for constructing alter¬ 
native basis functions. As noted in m, the main outstanding questions of PUM 
include the choice of the basis functions, the imposition of essential boundary 
conditions, and efficient numerical integration. The AES-FEM proposes a set 
of generalized Lagrange polynomial basis functions that also satisfy the parti¬ 
tion of unity. Therefore, AES-FEM effectively addresses these open problems 
in PUM. 
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Besides the aforementioned generalizations of FEM, it is also worth noting 
the least-squares finite element method (LSFEM) [11] . LSFEM uses the concept 
of least squares globally to minimize a global error. In contrast, AES-FEM uses 
least squares in a local sense for constructing basis functions. To use LSFEM, 
any higher order PDE must be decomposed into a system of first order PDEs 
first, so it does not preserve the framework of the standard FEM. 

Finally, we note the recent development of isogeometric analysis (IGA) |24j . 
which uses NURBS (Non-Uniform Rational B-Splines) or T-splines as basis func¬ 
tions instead of the standard FEM basis functions. These methods can deliver 
high accuracy over very coarse meshes and can be advantageous for problems 
that can benefit from high-degree continuity, such as thin-shell modeling. How¬ 
ever, IGA does not alleviate the dependency on mesh quality, since NURBS in 
effect impose stronger requirement on mesh quality than the standard FEM. 

2.3 Generalized Finite Difference and Weighted Least Squares 

The finite difference methods and finite element methods are closely related to 
each other. On structured meshes, the equivalence of these methods can be es¬ 
tablished in some special cases. While finite element methods were developed to 
support unstructured meshes from the beginning, the finite difference methods 
can also be generalized to unstructured meshes. These generalizations are often 
collectively referred to as generalized finite difference (or GFD) methods. The 
earlier GFD methods were based on polynomial interpolation; see e.g., na, na, 
and |26| . These methods construct a local multivariate polynomial interpola¬ 
tion by requiring the number of points in the stencil to be equal to the number 
of coefficients in the interpolant. However, due to the irregular distribution of 
points, the resulting Vandermonde matrices are often singular or ill-conditioned. 

More general than an interpolant are the least squares or weighted least 
squares approximations, which allow more points in the stencil than the number 
of coefficients of a polynomial. Some earlier examples of least-squares approx¬ 
imations include m and na, which attempted to improve the conditioning 
of interpolation-based GFD. More recently, the least-squares-based GFD have 
been analyzed more systematically by Benito et al. US! [571, and it been suc¬ 
cessfully applied to the solutions of parabolic and hyperbolic PDEs [UlllHlISZj 
and of advection-diffusion equation |29| . It has also been utilized in the weak- 
form of the Poisson equation, under the name meshless finite difference method 
(T6l[30|. In these methods, a weighting scheme is often used to alter the norm 
that is being minimized, and hence the name weighted least squares (WLS). 
These weights serve a different role from those in the moving least squares [20] . 
DEM, and EFGM. They do not need to have a compact support, and do not 
even to be defined by a continuous function in general. 

Even though the least-squares approximations tend to be better conditioned 
than their interpolation-based counterparts, ill-conditioning may still occur for 
a given set of points. To overcome the issue. Jiao et al. utilized adaptive stencils 
coupled with column and row scaling, QR factorization with column pivoting, 
and a condition number estimator, which effectively guarantee the conditioning 
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and hence the stability of approximations based on WLS [SllSa- InAES-FEM, 
we extend this previous work to construct the generalized Lagrange polynomial 
basis functions, and then utilize these basis functions in the weak form of the 
finite element methods. 


3 Formulation of AES-FEM 

The main idea of AES-EEM is the use of an alternative set of basis functions. 
For the basis functions, we propose the use of a set of generalized Lagrange 
polynomial basis functions (GLPBF) computed using a weighted least squares 
formulation. We use the standard FEM (hat) basis functions for the weight 
functions (a.k.a. test functions). In this section, we describe the weighted resid¬ 
ual formulation of FEM, define generalized Lagrange polynomial basis functions 
based on weighted least squares, and then describe AES-FEM. 

3.1 Weighted Residual Formulation of FEM 

We will approach the formulation of the finite element method through the lens 
of the weighted residual method for solving PDFs. The main idea of the for¬ 
mulation is to consider the unknown solution as a linear combination of basis 
(trial) functions and then select the coefficients such that the residual is orthog¬ 
onal to a set of weight (test) functions. Depending on the choice of the weight 
functions, one will derive different numerical methods, such as the collocation 
method, the least squares method, and the Galerkin method. Details about 
weighted residual and FEM can be found in [3311311133] • In the following, we 
give a brief overview of weighted residuals for completeness. 

Consider a linear differential operator C defined on a bounded, simply- 
connected domain O, with outward unit normal vector n. Denote the boundary 
ofDasP = ri:iUrAf, where r_D and Pat are disjoint sets on which Dirichlet 
and Neumann boundary conditions are specified, respectively. We want to hnd 
a function u such that 

Cu = f (I) 

subject to the boundary conditions 

Ou 

u = g on Td and 7 — = h on Pat. ( 2 ) 

on 

Eq. Ilj is the strong form of the PDE. In the weighted residual formulation, 
we use the weak form based on a set of weight functions dr = {ipi,... ,ipn}, by 
requiring the residual Lu — / to be orthogonal to ipi, i.e., 

[ {Cu -f)dV = 0. (3) 

Jq 

To approximate u, let $ = {4>i, be a set of basis functions, and we 
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define an approximation 


n 



( 4 ) 


Substituting (|3]) into the weak form ([3]) and rearranging the equations, we then 
obtain 





(5) 


At this point for simplicity, let us consider the Poisson equation with Dirichlet 
boundary conditions, for which the weak form is given by 



( 6 ) 


Substituting dU into m, we obtain 





(7) 


The finite element method uses integration by parts to reduce the order of 
derivatives required by ([7]). If tpi has weak derivatives and satisfies the con¬ 
dition '0i|r_D = then after integrating by parts and imposing the boundary 
conditions, we arrive at 



( 8 ) 


Taking ([5]) over the n weight functions, we obtain the linear system 


(9) 


Ku = g, 


where K is the stiffness matrix and g is the load vector, with 




( 10 ) 


and 


If the weight functions are chosen to be the same as the basis functions, then 
we will arrive at the Galerkin method. In this paper, we introduce a new set of 
basis functions based on weighted least squares and we use the standard linear 
FEM “hat functions” as the weight functions. 

3.2 Weighted Least Squares Approximations 

In this subsection, we review numerical differentiation based on weighted least 
squares approximations, as described in [siins- Similar to interpolation-based 
approximations, this method is based on Taylor series expansion. Let us take 
2D as an example, and suppose f{u) is a bivariate function with at least d + 1 


7 


continuous derivatives in some neighborhood of Uq = (0,0). Denote cji- = 
g^jgy/o Then for any u in the neighborhood, / may be approximated to 

the (d + l)st order accuracy about the origin Uq as 


d j+k=p j 

/N = E E (11) 

p^O j,k>0 ■ TT 

V. ^ remainder 

Taylor Polynomial 

Analogous formulae exist in ID and 3D. The derivatives of the Taylor polynomial 
are the same as those of / at Uq up to degree d. Therefore, once we have 
calculated the coefficients for the Taylor polynomial, finding the derivatives of 
/ at Uq is trivial. We proceed with calculating the coefficients as follows. 

For any point Uq, we select a stencil of m nodes from the neighborhood 
around Uq . Stencil selection is described further in Section [S] We do a local pa¬ 
rameterization of the neighborhood so that Uq is located at the origin (0,0) and 
the coordinates of the other points are given relative to Uq. Then substituting 
these points into d), we obtain a set of approximate equations 


d j+k=p 

^ 


p=0 j,k>0 




j\k\ 


fi, 


( 12 ) 


where fi = f {ui) and the Cjk denote the unknowns, resulting in an mxn system. 
There are n = {d+ \){d + 2)/2 unknowns in 2D and n = {d+ l){d + 2){d + 3)/6 
unknowns in 3D. Let V denote the generalized Vandermonde matrix, c denote 
the vector of unknowns (i.e., the Cjk) and f denote the vector of function values. 
Then we arrive at the rectangular system 

Vc^f. (13) 


Let us now introduce some notation to allow us to write the Taylor series 
in matrix notation before we proceed with our discussion of solving (I13II . Let 
T’lf\x) denote the set of all /c-dimensional monomials of degree d and lower, 
stored in ascending order as a column vector. If no ambiguities will arise, we 
will use V in place of For example, for second degree in 2D we have 

(x) = [l X y xy (14) 

Let D be a diagonal matrix consisting of the fractional factorial part of the 
coefficients, i.e. in (fTTl) . For example, for second degree in 2D we have 

D = diag ^1, 1, 1, 1, . (15) 

Then we may write the Taylor series as 

f[x) = c^DV {x ). 
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(16) 







To solve m, we use a weighted linear least squares formulation [33], that 
is, we will minimize a weighted norm (or semi-norm) 


imn\\Vc- f\\^^ = mm\\W {Vc- (17) 


where W is a,n m x m diagonal weighting matrix. The entries of W assign 
weights to the rows of matrix V. Specifically, if we denote the diagonal entries 
of W as Wi, then row i is assigned weight Wi. These weights can be used to 
prioritize the points in the system: we assign heavier weights to the nodes that 
are closer to the center point. By setting a weight to zero (or very close to zero), 
we may also filter out outliers or other undesirable points. Note that for a given 
node, the weighting matrix W is constant. 

If / is in the column space of V, then the solution of the linear system is 
not affected by a nonsingular weighting matrix. However, if / is not in the 
column space, which is often the case, then different weighting schemes can lead 
to different solutions. Choosing a good weighting matrix is application specific. 
For quadratic approximations, we compute the weights as follows. Let h denote 
the maximum radius of the neighborhood, that is 


Then 


Wi = 


max 

l<2<m 

(18) 

("-ik+d . 

(19) 


where e is a small number, such as e = 0.01, for avoiding division by zero. 

After the weighting matrix has been applied, we can denote the new system 
as 

Me ft! /, where M = WV and / = Wf. (20) 


This resulting system may be rank-deficient or ill-conditioned. This is a chal¬ 
lenge that GFD researchers have been dealing with since the 1970s [T3|. The 
ill-conditioning may arise from a number of issues including poor scaling, an 
insufficient number of nodes in the neighborhood, or a degenerate arrangement 
of points. We resolve these issues with neighborhood selection, discussed in 
Section O We can address the scaling issue with the use of a diagonal scaling 
matrix S. Let aj denote the jth column of an arbitrary matrix A. A typical 
choice for the jth entry of S is either 1/ ||aj||2, which approximately minimizes 
the 2-norm condition number of AS |36|, or 1/ 11% [37]. Using exact arith¬ 
metic, the matrix S does not affect the solution, but it can significantly improve 
the conditioning and thus the accuracy in the presence of rounding errors. After 
applying the scaling matrix to WV, the problem becomes 


min 

d 


Vd-f 


where V = WVS = MS and d = S-^c. 


( 21 ) 


Conceptually, the solution to the above problem may be reached through the 
use of a pseudoinverse. We will have 

d=V^f where = (v^V^ ^ ■ (22) 
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However, since the resulting system may still be rank-deficient or ill-conditioned, 
we solve it using QR factorization with column pivoting, as discussed in Sub¬ 
section 15.41 Finally, we get the vector of partial derivatives for the Taylor 
polynomial 

c = Sd. (23) 

3.3 Generalized Lagrange Polynomial Basis Functions 

We now define basis functions based on weighted least squares. Note that the 
standard finite element methods use piecewise Lagrange polynomial shape func¬ 
tions, which have two especially important properties: the coefficients of the 
basis functions have the physical meaning of the function values or their ap¬ 
proximations at the nodes, and the basis functions form a partition of unity. 
We refer to the two properties as function value as coefficient and partition of 
unity, respectively. These properties are desirable in ensuring the consistency 
of the method based on these basis functions and also for the ease of imposing 
Dirichlet boundary conditions. However, the traditional concept of the Lagrange 
basis functions is interpolatory, so they are not applicable to least squares. We 
now generalize this concept, so that it can be applicable to least-squares-based 
basis functions. 

Definition 1. Given a set of degree-d polynomial basis functions {4>i}, we say 
it is a set of degree-d generalized Lagrange polynomial (GLP) basis functions if: 

1. / (xi) (fi approximates a function f to O in a neighborhood of 

the stencil, where h is some characteristic length measure, and 

2 . = 

We now define a set of GLP basis functions based on weighted least squares. 

Given a stencil {xt}, we follow the procedure in Subsection 13.21 When 
computing the jth basis function (fj, let f = ej, where Sj is the jth column of 
the identity matrix. Following and (1^5)) . we have 

c = SV^Wej. (24) 

Thus for the jth basis function, the vector c is exactly the jth column of 
SV W. We define a set of basis functions as 

§> = (^SV^wY DV. (25) 

For a more concrete example, if we denote the elements of as we can see 
that the jth basis function for degree 2 in 2D is 

(fj = Wj ^aijSi -I- a2jS2X + asjSsy + 04 ^ 34 ^^ + a^jS^xy + . (26) 

Note that Wj is a constant scalar, so cfj is a polynomial. The basis functions in 
(I 25 I) are an example of GLP basis functions. We summarize this key feature in 
the following theorem. 
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Theorem 2. The basis functions based on weighted least squares as defined in 
are generalized Lagrange polynomial basis functions. 


We shall postpone the proof of this theorem to Section SI where we will also 


analyze the accuracy and stability of finite element discretization based on these 
basis functions. In the following, we will finish the description of AES-FEM. 


3.4 Description of AES-FEM 


Starting with the weighted residual formulation for FEM from ([5]), we propose 
using GLP basis functions for the basis functions in the weak form and using the 
traditional hat functions for the weight functions. More specifically, for a given 
node i and its associated weight function ipi , a specific set of GLP basis functions 
{(t>j} is constructed from a weighted stencil {Xi,Wi) of n neighboring vertices 
centered at node i. This weight function '0^ and its associated set of GLP basis 
functions {(j)j} are used to compute the ith row of the stiffness matrix, as given 
by cni). Note that, because a different set of basis functions is associated with 
each weight function, the basis functions on a given element differ row-to-row 
in the stiffness matrix. 

For AES-FEM, when using the 1-ring neighborhood of a vertex as the sten¬ 
cil, we can use quadratic GLP basis functions, which have the advantage of 
decreased dependence on element quality and improved accuracy over standard 
finite element with linear basis functions, while virtually preserving the sparsity 
pattern of the stiffness matrix, as we will demonstrate in Section [6l 

In terms of the computation of the load vector, we can use either the FEM 
or AES-FEM basis functions. We refer to the AES-FEM with these two options 
as “AES-FEM 1” and “AES-FEM 2,” respectively. Additional implementation 
details will be given in Subsection 15.41 


4 Analysis of AES-FEM 


In this section, we analyze the consistency and stability of AES-FEM. We start 
by explaining the applicability of Green’s identity to the GLP basis functions. 
We will then prove that basis functions in (1^51) are GLP basis functions and 
discuss the consistency and stability of AES-FEM. 

4.1 Green’s Identity and Integration by Parts 

For a given weight function 0, the GLP basis functions are continuously differ¬ 
entiable everywhere in within the domain of integration. Hence, in a variational 
sense, the GLP basis functions still allow one to use Green’s identities to formu¬ 
late a weak form. Let (j)j be any GLP basis function computed from a weighted 
stencil, and ipi be a classical FEM shape function with compact support con¬ 
tained within the set O. Therefore, for any partial derivative operation d, it 
follows that 



( 27 ) 
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Because of this property, the weak-form formulation with the GLP basis func¬ 
tions in AES-FEM is mathematically sound. 

In addition, since the finite-element basis functions tpi vanish along the 
boundary, after integration by parts, we have the identities 

[ dV = - [ VV-* ■ dV. (28) 

Jq. Jq 

Therefore, we can reduce the order of derivatives similar to the classical FEM, 
without introducing additional boundary integrals to the computation. 

4.2 Generalized Lagrange Polynomial Basis Functions 

We now show that the basis functions in (1251) are indeed generalized Lagrange 
polynomial basis functions, which follows from the two lemmas below. 

Lemma 3. Let {xi} be a stencil with m nodes and stencil diameter h. Let {4>j} 
be the complete set of basis functions of degree up to d as defined by i25\) on this 
stencil. Given an arbitrary function f, define the following approximation of f 

m 

f'^{x) = J2fjM^) = f^^^ ( 29 ) 

i=i 

where fj = f {xj) and f = [fi f 2 ■■■ fm]^- If the rescaled matrix V has a 
bounded condition number, then f^ approximates f to O . In addition, 

given any degree-k differential operator D, if f is continuously differentiable up 
to degree k, then Vf^ approximates T>f to O 

Proof. First, we show that the approximation is equivalent to directly solving 
a weighted least squares problem for a local polynomial fitting of /. Using the 
method described in Subsection l3.2l we get the coefficients c = SV ^Wg where 
9 =[fif 2 • ■ • /m]^- Thus, the local polynomial fitting of / is 

/ (a;) « c^DV 

= (^SV^WgY DV 

= (SV'^W'Y DV 

= g^^ 

= f (*)■ 

It follows from Theorem 4 in [21] that approximates f to O and Vf^ 

approximates Vf to O for any degree-fc differential operator V. □ 

Lemma 4. The basis functions in i25\) form a partition of unity, i. e., 4^j (^) 

1 . 
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Proof. We will show that on a given stencil {xi}, the GLP basis functions {(j)j} 
form a partition of unity. Let V be the generalized Vandermonde matrix for 
the given stencil. For example, for second order expansion in 2D, we have 


1 xi yi ... ij/f 

1 X2 2/2 • ■ • Ivj 

1 Xm ym • ■ • 2 ^ 7 ^ 


(30) 


For a given function /, using the truncated Taylor series, we have 

Vc = g, (31) 

where c is the vector of partial derivative values. Applying the diagonal row 
weighting matrix W and the diagonal column scaling matrix S, we have 

V (S'-^c) = Wg where V = WVS. (32) 

This is a least squares problem, and the solution for c may be reached through 
the use of a pseudoinverse, 

c = SV^Wg. (33) 

For the jth GLP basis function, we have g^ = [0... 1... 0]^, where the 1 is in 

the jth position, and hence the columns of SV^W multiplied by the Taylor 
constants D give the coefficients for the GLP basis functions. This implies that 
the entries of the ith row of SV W correspond to the coefficients of the ith 
terms in the set of basis functions. 

To finish the proof, it suffices to show that the sum of the entries in the first 
row of SV W is 1, and the sum of the entries in any other row is 0. Let vector 
w be the diagonal entries of W. Every entry of the first column of V is equal 

to 1, thus the first column of V is then siw. Denote the ith row of V as .j. 

The sum of the entries of the ith row of SV W is Since V is a left 

inverse of V, we have V~^V = I, and hence 


~+T 


1=1 
2 < i < n 


(34) 


Therefore, the GLP basis functions form a partition of unity. □ 

From the above lemmas, the basis functions in (l25)l satisfy both the prop¬ 
erties of function value as coefficient and partition of unity, and hence they are 
GLP basis functions, as claimed in Theorem [5J 
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4.3 Consistency of AES-FEM 

The accuracy of the AES-FEM depends on its consistency and stability. We 
hrst consider the consistency of AES-FEM, in terms of its truncation errors in 
approximating the weak form ([3]) by ([5]) using the GLP basis functions. In a 
nutshell, the consistency of the AES-FEM follows directly from Lemma [31 For 
completeness, we consider a specific example of solving the Poisson equation. 
The analysis for other PDFs can be derived in a similar fashion. 

Theorem 5. Suppose u is smooth and thus ||Vu|| is bounded. Then, when 
solving the Poisson equation using AES-FEM with degree-d GLP basis functions 
in for each ifi the weak form is approximated to 0{h‘^), where h is some 
characteristic length measure of the mesh. 

Proof. Let u be the exact solution on a mesh with mesh size h, and let 

n 

u = 

j=i 


denote the approximation to u using degree-d GLP basis functions. Because the 
test functions xpi vanishes along boundary, the weak form ([31) can be rewritten 
as 


- [ Vu-VipidV = I ififdV. 

Jq Jq 

When using degree-d GLP basis functions, it follows from Lemma [3] that 


||Vu-V m|| =0{h‘^) 


within each element. Under the assumption that u is twice differentiable, Vm is 
bounded, and hence 


/ (Vm - Vm) • dV 

JQ. 


= O (h^) . 


Furthermore, if / is approximated by degree-d GLP basis functions as 


i=i 


then we have 




dV 


O . 


□ 


More specifically, when using quadratic GLP basis functions, the truncation 
errors are second order in the stiffness matrix. The truncation errors in the 
load vector is third order for AES-FEM 2 when / is also approximated using 
quadratic GLP basis functions, but it is second order for AES-FEM 1 when / is 


14 






approximated using linear FEM hat functions. Like most other PDE methods, 
as long as the method is stable, the rounding errors do not dominate the trun¬ 
cation errors, and there is no systematic cancelation of truncation errors, we 
expect the solution to converge at the same rate as the local truncation errors, 
as we will demonstrate numerically in Section [51 

4.4 Stability 

Eor elliptic PDEs, the stability of a method depends on the condition number of 
its coefficient matrix, which can affect the performance of iterative solvers and 
the accuracy of the solution. It is well known that the traditional finite element 
method may be unstable for poorly shaped meshes [3B], and some meshless 
methods may also suffer from instability when two points nearly coincide. AES- 
EEM avoids these potential instability issues. 

As a concrete example, let us consider the Poisson equations with Dirichlet 
boundary conditions, whose coefficient matrix is the stiffness matrix. It is well 
known that the condition number of the stiffness matrix is proportional to h~^, 
where h is some characteristic length of the mesh [39]. However, if the condition 
number is significantly larger, then the method is said to be unstable, which can 
happen due to various reasons. 

The ill-conditioning of any local stiffness matrix may lead to poor scaling 
and in turn ill-conditioning of the global stiffness matrix. This is owing to the 
following fact, which is given as Theorem 2.2.26 in BOj. 

Proposition 6. For any matrix A G m > n, its condition number 

in any p-norm, denoted by Hp (A ), is hounded by the ratio of the largest and 
smallest column vectors in p-norm, i.e., 


(^) — 


maxi<i<„ llajp 
mini<j<n llajllp’ 


where denotes the kth column vector of A. 


(35) 


The above fact offers an intuitive explanation of a source of ill-conditioning in 
traditional finite element methods due to poorly shaped elements: poorly shaped 
elements may lead to unbounded large entries in local stiffness matrices, so the 
column norms of the global stiffness matrix would vary substantially, and in turn 
the global stiffness matrix is necessarily ill-conditioned. In the context of AES- 
EEM, there can be two potential sources of local instability due to poor scaling. 
Eirst, the unnormalized local Vandermonde system given in (1131) is in general 
very poorly scaled. We resolved this by normalizing the Vandermonde system 
to avoid poor scaling. Second, the normalized Vandermonde system may still be 
ill-conditioned occasionally, when a stencil is degenerate or nearly degenerate, 
which could lead to unbounded large values in the local stiffness matrix. We 
resolved this issue by using QR with column pivoting and condition-number 
estimation. 
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Even if the local stiffness matrices are bounded, the global stiffness matrix 
may still be ill-conditioned due to linearly dependent rows or columns. This is 
the potential source of instability for some meshless methods when two points 
nearly coincide; the two points may share the same stencil and basis functions, 
so that the rows or columns corresponding to the two points would be nearly 
identical. Therefore, these meshless methods also require good point distribu¬ 
tions. In AES-FEM, we utilize the mesh topology to construct the stencil, as 
we will describe in Section [5l This ensures that no two vertices share the same 
stencil unless there are coincident points, and hence it gives a strong guarantee 
that the rows in the global stiffness matrix are linearly independent. 

The aforementioned reasons are the most common causes of instability for 
solving elliptic PDEs. Another source of instability is a cluster of coincident 
points or inverted elements, which rarely happen in practice, and hence we 
defer their treatments to future work. As we will demonstrate numerically in 
Section [SI by resolving these instabilities, AES-FEM produces well-conditioned 
stiffness matrices for meshes, even with very bad quality elements or point 
distributions. 


5 Implementation 

We discuss the practical aspects of the implementation of AES-FEM in this 
section. We start with a discussion of the utilized mesh data structure and then 
explain how this enables quick and efficient neighborhood selection. Finally, the 
algorithms are presented and runtime is analyzed. 

5.1 Data Structure 

We use an Array-based Half-Facet (AHF) data structure [IT] to store the mesh 
information. In a d-dimensional mesh, the term facet refers to the (d — 1)- 
dimensional mesh entities; that is, in 2D the facets are the edges, and in 3D the 
facets are the faces. The basis for the half-facet data structure is the idea that 
every facet in a manifold mesh is made of two half-facets oriented in opposite 
directions. We refer to these two half-facets as sibling half-facets. Half-facets 
on the boundary of the domain have no siblings. In 2D and 3D, the half-facets 
are half-edges and half-faces, respectively. We identify each half-facet by a two 
tuple: the element ID and a local facet ID within the element. In 2D, we store 
the element connectivity, sibling half-edges, and a mapping from each node to an 
incident half-edge. In 3D, we store the element connectivity, sibling half-faces, 
and a mapping from each node to an incident half-face. For an example of a 
2D mesh and the associated data structure, see Figure |TJ This data structure 
allows us to do neighborhood queries for a node in constant time (provided the 
valance is bounded). For additional information about the AHF data structure, 
see [TT] . 
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Figure 1: An example of half edges and associated data structure. 

5.2 Neighborhood Selection 

The use of the AHF data structure allows us to quickly find the neighborhood 
of a node. We use the concept of rings to control the size of the neighborhood. 
The 1-ring neighbor elements of a node are defined to be the elements incident 
on the node. The 1-ring neighborhood of a node contains the nodes of its 1- 
ring neighbor elements m- Most of the time, when using GFD with second 
order basis functions or when constructing second order GLP basis functions, 
the 1-ring neighborhood of a node supplies the appropriate number of nodes. 
If the valance is low, it might be necessary to further expand and collect more 
nodes for the neighborhood. Therefore, for any integer fc > 1, we define the 
{k -\- l)-ring neighborhood as the nodes in the /c-ring neighborhood plus their 
1-ring neighborhoods. 

As k increases, the average size of the fc-ring neighborhood grows very 
quickly. The granularity can be fine-tuned by using fractional rings. In 2D 
we use half rings, which are defined in (31]; for any integer k > 1 the {k -\- 1 / 2 )- 
ring neighborhood is the fc-ring neighborhood plus the nodes of all the faces that 
share an edge with the fc-ring neighborhood. See Figure [5] for a visualization of 
rings and half-rings in 2D. We extend this definition to 3D and introduce 1/3- 
and 2 / 3 -rings. For any integer fc > 1, the (fc -I- ^/^Yring neighborhood contains 
the fc-ring neighborhood plus the nodes of all elements that share a face with 
the fc-ring neighborhood. The (fc -I- ‘^Y)-ring neighborhood contains the fc-ring 
neighborhood plus the nodes of all faces that share an edge with the fc-ring 
neighborhood. 

Note that for 2D triangular and 3D tetrahedral meshes, the 1-ring neighbor¬ 
hood typically has enough points for constructing quadratic GLP basis func¬ 
tions. Therefore, the stiffness matrix from AES-FEM has a similar sparsity 
pattern to that from standard FEM with linear shape functions. However, 
when the 1-ring neighborhood is too small, the extended stencil with a larger 
ring allows AES-FEM to overcome mesh-quality dependence and improve its 
local stability. 
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Figure 2: Examples of 2D stencils with 1-ring, 1.5-ring, 2-ring, and 2.5-ring 
neighborhoods of center node (in solid black). 

5.3 Stable Computation of GLP Basis Functions 

Even with proper neighborhood selection, it cannot be guaranteed that the 
least-squares problem in (l?T]) will be well-conditioned. Hence, it is critical to 
use a robust method that ensures the accuracy and stability of the approximate 
solutions. 

Note that one standard technique in linear algebra for solving rank-deficient 
least-squares problems is truncated singular value decomposition (TSVD) [36]. 
The TSVD is not recommended here, because it can result in the loss of partition 
of unity of the basis functions. The reason is as follows. When truncating SVD, 
one truncates any singular value aj that is smaller than ecri, where cri is the 
largest singular value and e is some small positive value, such as 10“^. These 
singular values may be necessary for computing the constant terms in the GLP 
basis functions, and hence their loss can result in a set of basis functions that 
lack the partition of unity and in turn may compromise convergence. 

We avoid the above issue by using truncated QR factorization with column 
pivoting (QRCP). When performing QRCP, one can find a numerical rank r 
of matrix R so that the condition number K{Ri-,r,i:r) < 1/e, where e is some 
small positive value, such as 10“'^. This is elaborated in the below discussion 
of Algorithm [T] If r is less than the size of R, then any diagonal entry of R 
in position r -|- 1 or greater is truncated, thus truncating the (r + l)th and 
subsequent columns of Q. In QRCP, we require the first column not to be 
permuted, and this ensures the resulting basis functions satisfy the property of 
partition of unity. 

In Algorithmlll we present the procedure for initializing the generalized Van¬ 
dermonde matrix V and factoring it using QRCP. The generalized Vandermonde 
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Algorithm 1 Initialization of a Generalized Vandermonde Matrix 

function: initiate_GVM 

input: 1 . Xk'. local coordinates of stencil 

2. w. vector of row weights 

3. p: desired degree for V 

4. e: tolerance for rank deficiency 

output: struct gvm: with fields TV, S, Q, R, P, r (estimated rank) 

1; create generalized Vandermonde matrix V from local coordinates Xk 
2: determine column scaling matrix S 
3: TV t— diag(u;) 

4: V ^ WVS 
5; solve VP ~ QR 

6; estimate rank r from R so that r = max{j|cond (.Riii.i:*) < 1/e} 


matrix is formed from the local coordinates of the stencils and is scaled by the 
column scaling matrix S and the row scaling matrix TV. The resulting matrix 
is then factored using QRGP. We use a variant of Householder triangularization 
|36| since this procedure is more efficient and stable than alternatives (such as 
Gram-Schmidt orthogonalization). When implementing this procedure, the QR 
factorization of V can overwrite V. The jth Householder reflection vector is 
of size n — j + 1. By requiring the first element of the vector to be positive, 
the first element may be reconstructed from the other elements, and thus only 
n — j entries are required to store the jth Householder reflection vector. The 
Householder vectors are stored in the lower triangular part of V and the R 
entries are stored in the upper part. The permutation matrix P is stored in a 
permutation vector. 

In addition to computing the QR factorization of the generalized Vander¬ 
monde matrix, an estimation of the numerical rank of R is also computed in 
Algorithm [T] The rank is important for ensuring the overall stability of other 
algorithms that use this initialization step. In order to estimate the numerical 
rank, we estimate the condition numbers of the leading principal sub-matrices 
of i?, Ri-.r.i-.r’, and find the largest r such that k(i2i:r,i:r) £ 1/e where k is 
the estimated condition number and e is some given drop-off tolerance depend¬ 
ing on the degree of polynomials. Note that since the matrix R is small, the 
condition numbers in different norms differ by only a small factor. Therefore 
for efficiency, we estimate the condition number of R in the 1-norm using the 
algorithm described in |42| . 

Once the generalized Vandermonde matrix has been initialized, it may be 
used to construct generalized finite differentiation operators from the weighted 
least squares approximations, as described in Algorithm [21 The input for this 
algorithm is the output from Algorithm [T] and a vector a. The vector a = 
VPix) contains the values for some specified derivative V of the monomial 
basis functions P at point x. For example, let T) he Then in 2D, we 
have a = 'DP{x,y) = [0 0 1 0 a; and in 3D, we have a = VPix^y, z) = 
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Algorithm 2 Approximating Vf at given point x from WLS 
function diff_WLS 

inpnt: 1. struct gvm: with fields W, S, Q, R, P, r (estimated rank) 

2. coefficients a = W{x) 

outpnt: weights d, so that g = 'Df{x) for g containing f(xk) at stencil 
points 

1: a <— {P. S~^a; 

2: d i — R\'r l-r^5 
3: d -(r- WQ. ^a] 


[0 0100x0 2y0 0]^. The algorithm returns a vector of weights d so that 
d^g = 'Df{x) for a vector g = [fi f 2 ■■ ■ fm] containing the values of the 
function at the stencil points. Note that for a GLP basis function, the returned 
weights are the values of the specified derivative at the points in the stencil. 

In terms of the computational cost, the step that dominates Algorithm [1] is 
the QRCP factorization, which takes O (2mn^ — fri^) flops where V is m x n 
|36| . Here, m is the number of points in the stencil and n is the number of 
terms in the Taylor series expansion (for second order expansion n = 6 in 2D 
and n = 10 in 3D). As long as the valance is bounded, that is m is bounded, this 
algorithm is executed in a constant time. Compared to calculations based on the 
standard finite-element basis functions, which are tabulated, the computation 
based on the GLP basis functions is more expensive. This leads to higher cost 
of AES-FEM in assembling the stiffness matrix and load vector, as we discuss 
next. However, this cost is a small constant per element, and AES-FEM can be 
more efficient overall by delivering higher accuracy, as we will demonstrate in 
Section [51 

5.4 Assembly of Stiffness Matrix and Load Vector 

Algorithm 13] presents a summary of the AES-FEM procedure for assembling the 
stiffness matrix and load vector for a PDE with Dirichlet boundary conditions. 
Unlike the standard FEM procedure, we build the stiffness matrix row by row, 
rather than element by element. This is because the most computationally 
expensive part of the procedure is to compute the derivatives for the set of basis 
functions on each stencil. A weight function is nonzero only on the neighborhood 
around its corresponding node. Since the weight functions correspond to the 
rows, we assemble the stiffness matrix row by row, ensuring that we will only 
need to compute the derivatives for each neighborhood once. 

When computing a row of the stiffness matrix, the first step is to obtain 
the stencil of node k. This step is performed by utilizing the data structure 
presented in Subsection 15.11 and the proper size of the stencil is ensured by 
choosing the ring sizes adaptively. Next, the local coordinates are calculated for 
the points in the stencils and the row weights are computed. Using Algorithm[Tl 
the QR factorization of the generalized Vandermonde matrix is computed for 
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the neighborhood. Then for each element that contains node k, we perform 
the integration of the weak form in a manner that is similar to standard FEM. 
The element Jacobian is computed and used to find the local coordinates of the 
quadrature points. The derivatives of the weight function are computed, that 
is V'0i, at the quadrature points of the current element. Recall that the weight 
functions are the standard hat functions. Next the derivatives of the basis 
functions, that is V(^j, are computed at the quadrature points of the current 
element. The basis functions are the GLP basis functions and thus Algorithmic] 
is used. The value of the integral on the current element is computed and either 
added to the stiffness matrix or subtracted from the load vector, depending 
on whether the basis function corresponds to a node with Dirichlet boundary 
conditions. 

When computing the load vector, typically the entries h = J 'ipif dV are 
computed using a quadrature rule. One may evaluate / at the quadrature points 
in two ways. The first way is to use the standard procedure in FEM, i.e., to use 
the FEM basis functions and the values of / at the nodes of the element. Let 
M be the vector of EEM shape functions evaluated at quadrature point Xk and 
Seiem ^>6 the vector of the function values at the nodes of the element. Then, 
we have 


/ [xk) = M ■ 9^13^. 


(36) 


Alternatively, we may use GLP basis functions to interpolate the values of / at 
the quadrature points. We can approximate an arbitrary function using the set 
of GLP basis functions. In matrix notation, we have 



(37) 


where is the vector of function values at the nodes in the stencil and the 
vector "P (xk) has been evaluated at the quadrature point Xk- As mentioned 
earlier, we refer to the variant of AES-EEM using the former method of calculat¬ 
ing the load vector as AES-FEM 1 and refer to the latter variant as AES-FEM 
2 . 

We use Gaussian quadrature to perform the integration within each element. 
Eor quadratic GLP basis functions in 2D, we use a 1-point rule for stiffness 
matrix, and a 3-point rule for the load vector. In 3D, we use a 1-point rule 
for the stiffness matrix and a 4-point rule for the load vector. These rules are 
exact because the basis functions and their derivatives are quadratic and linear, 
respectively. 

It is worth noting that because of the properties of generalized Lagrange 
polynomial basis functions, Dirichlet boundary conditions may be imposed in 
AES-FEM in the same manner as in standard FEM. One does not need to use 
Lagrange multipliers or a penalty method. Additionally, the standard method 
for imposing Neumann boundary conditions may be used in AES-EEM. 

All the steps inside of the primary for-loop are executed in constant time, 
assuming that the size of each neighborhood is bounded. Therefore, the assem¬ 
bly of the stiffness matrix in AES-FEM has an asymptotic runtime of 0{n), 
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where n is the number of nodes in the mesh. When using AES-FEM 2, that 
is when WLS approximation is used to compute the approximation of / at the 
quadrature points, Algorithm is called again. While this function is constant 
in runtime, it has a large coefficient and thus takes longer than approximating 
the values of / using FEM (hat) basis functions. Therefore, the assembly time 
for AES-FEM 2 is longer than that for AES-FEM 1, as can be seen in Section [51 

Once the stiffness matrix and the load vector are assembled, we use the 
generalized minimal residual method (GMRES) [13] with a preconditioner to 
solve the linear system. GMRES is a standard Krylov subspace method for 
iteratively solving a sparse, nonsymmetric linear system. Specifically, we use the 
MATLAB implementation of GMRES. Another option for solving this system 
would be multigrid methods |44| . and we will explore their use in future work. 
As a preconditioner, we use incomplete LU factorization (ILU) in 2D and Gauss- 
Seidel in 3D. 

Finally, for completeness, we include Algorithm H] to summarize the GFD 
procedure for solving PDFs with Dirichlet boundary conditions. We use this 
algorithm primarily for comparison with the AES-FEM algorithm. We can see 
that for GFD, we need to use Algorithm [2] to compute the weights once for 
each non-Dirichlet node in the mesh; that is, if there are n non-Dirichlet nodes. 
Algorithm is called n times. For AES-FEM, for a given node, we use this 
algorithm once for every element containing that node. Thus if every node has 
a neighborhood of k elements and there are n non-Dirichlet nodes, we call the 
algorithm kn times. Therefore, the assembly time is longer for AES-FEM than 
for GFD, as we will see in Section ID 

6 Numerical Results 

In this section, we compare the accuracy, efficiency, and element shape quality 
dependence of AES-FEM 1, AES-FEM 2, FEM with linear basis functions, and 
GFD with quadratic basis functions. We compare these four method because the 
sparsity pattern of the coefficient matrix is nearly identical for all the methods. 
The sparsity pattern determines the amount of storage necessary and also the 
computational cost of vector-matrix multiplication. The errors are calculated 
using the discrete L 2 and Loo norms. Let u denote the exact solution and let u 
denote the numerical solution. Then, we calculate the norms as 

L 2 (error) =^y lu — and Loo (error) = max lu — u|. (38) 

For each series of meshes of different grid resolution, we calculate the average 
convergence rate as 

/ error of mesh 1\ /, / d I nodes in mesh 1 \ 

convergence rate = - log 2 - - -^ / log 2 \ ——^^ , 

\ error of mesh 4 / / \ V nodes m mesh 4 j 

(39) 

where d is the spacial dimension. 
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Algorithm 3 Building a Stiffness Matrix and Load Vector using AES-FEM 
function: aes_feni 

inpnt: 1. x, elem, opphfs, vh2f: mesh information 

2. p: desired degree for GLP functions 

3. e: tolerance for rank deficiency 

4. AESFEMl: boolean for AES-EEM 1 or AES-FEM 2 

5. isDBC: flags for Dirichlet boundary conditions 
outpnt: stiffness matrix K and load vector b 

1: for each node without Dirichlet boundary conditions do 
2 : obtain neighborhood of node 

3; calculate local parameterization Xk and row weights w for neighborhood 
4; aes_gvm ^ initiate_GVM(a;fc, w, p, e) 

5; obtain local element neighborhood 
6: for each element in local neighborhood do 

7: calculate element Jacobian and local coordinates of quad-points 

8; calculate derivatives of FEM shape functions at quad-points 
9; a ^ VPlx) where 'D'P{x) is defined by the PDE we are solving 
10: GLPderivs ^ diff_WLS(aes_gvm, a) 

11: for each node in neighborhood do 

12: if not Dirichlet BG node then 

13: add integral to appropriate stiffness matrix entry 

14: else 

15: subtract integral from load vector 

16: end if 

17: end for 

18: if AESFEMl then 

19: calculate load vector over current element using FEM approximations 

for quad-points 
20: else 

21: calculate load vector over current element using GLP approximations 

for quad-points 
22: end if 

23: end for 

24: end for 
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Algorithm 4 Constructing a GFD coefficient matrix 

function: gfd 

input: 1. X, elem, opphfs, vh2f: mesh information 

2. p: desired degree for GFD functions 

3. isDBC: flags for Dirichlet boundary conditions 
output: GFD matrix K and vector b 

1: for each node without Dirichlet boundary conditions do 
2; obtain neighborhood of node 

3; calculate local parameterization and row weights for neighborhood 
4; gfd_cvm ■«— initiate_GVM(a:fc, w, p, e) 

5: a e- W{x) where 'D‘P{x) is defined by the PDF we are solving 

6: GFDderivs t— diff_WLS(gfd_cvm, a) 

7: for each node in local neighborhood do 

8; if not BG node then 

9; enter value in matrix 

10: else 

11: subtract from RHS vector 

12: end if 

13: end for 

14: end for 


6.1 2D Results 

In this section, we present the results of our 2-dimensional experiments. We 
use two different series of meshes, each series with 4 meshes. The first series 
of meshes (referred to collectively as “mesh series 1”) is generated by placing 
nodes on a regular grid and then using MATLAB’s Delaunay triangularization 
function to create the elements. The meshes range in size from 64 x 64 to 
512 X 512 nodes. On the most refined mesh, the minimum angle is 45 degrees 
and the maximum angle is 90 degrees. The maximum aspect ratio is 1.41, where 
a triangle’s aspect ratio is defined as the ratio of the length of longest edge to the 
length of the smallest edge. The second series of meshes (referred collectively 
as “mesh series 2”) is generated by using Triangle [45]. The number of nodes 
for each level of refinement is 4,103, 16,401, 65,655, and 262,597, respectively, 
approximately the same as those in series 1. On the most refined mesh, the 
maximum angle is 129.6 degrees and the minimum angle is 22.4 degrees. The 
maximum aspect ratio is 2.61. See Figure [3] for a visualization of the types of 
meshes used. 
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Figure 3: The mesh on the left is representative of the meshes used in series 1. 
The mesh on the right is representative of the meshes used in series 2. Note that 
the meshes above are coarser than the meshes used in computations so that the 
details can be seen clearly. 


6.1.1 Poisson Equation 

The first set of results we present is for the Poisson equation with Dirichlet 


boundary conditions on the unit square. That is, 

-V^u = f inf2 = [0,l]2, (40) 

u = g on dVL. (41) 

We consider the following three analytic solutions: 

ui = 16a;(l - a:)y(l - y), (42) 

U 2 = cos(7ra:) cos(7ry), (43) 

U 3 = 5 — sinh(7ra;) cosh(7ry). (44) 

smh TT cosh tt 


The Dirichlet boundary conditions are obtained from the given analytic solu¬ 
tions. The boundary conditions for ui are homogeneous and they are non- 
homogenous for U 2 and U 3 . 

The Too and L 2 norm errors for ui on mesh series 1 are displayed in Figure 
01 One can see that the two graphs are fairly similar; this is true for U 2 and U 3 
as well and thus we show only the Too norm errors for these two problems; see 
FigurejS] GFD is the most accurate for ui and U 2 and AES-FEM 2 is the most 
accurate for U 3 . EEM is the least accurate in all three cases. 

Eor mesh series 2, the Too and T 2 norm errors for ui can be seen in Figure 
inland the Too norm errors for U 2 and U 3 can be seen in Figure [71 On this mesh 
series, AES-FEM 2 has the lowest error for ui and U 2 - For U 3 , the errors for 
GFD and AES-FEM 2 are very similar. 
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Figure 4: The errors for 2D Poisson equation on mesh 1 for ui. The errors were 
computed using the norm (left) and the L 2 norm (right). 




Figure 5: The Loo norm errors for the 2D Poisson equation on mesh 1 for U 2 
(left) and M3 (right). 




Figure 6: The errors for 2D Poisson equation on mesh 2 for ui. The errors were 
computed using the Loo norm (left) and the L2 norm (right). 
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Figure 7: The Loo norm errors for the 2D Poisson equation on mesh 2 for U 2 
(left) and M 3 (right). 


6.1.2 Convection-Diffusion Equation 

We consider the convection-diffusion equation with Dirichlet boundary condi¬ 
tions on the unit square. That is, 

-V^u -h c • Vm = / in D, (45) 

u = g on on. (46) 

We take c = [1,1]^ for all of our tests and we consider the same analytic solu¬ 
tions as for the Poisson equation. Again the boundary conditions are obtained 
from the given analytic solutions. 

The Loo and L 2 norm errors obtained from the convection-diffusion equation 
on mesh series 1 with ui are presented in Figure [51 The Loo norm errors for 
U 2 and U 3 on mesh series 1 are in Figure |9l For all three problems, AES-FEM 

and GFD are both more accurate than linear FEM. For ui, GFD is the most 

accurate. For U 2 , the most accurate method is either AES-FEM 2 or GFD 
depending on the level of refinement. For U 3 , AES-FEM 2 is the most accurate. 

On mesh series 2, AES-FEM 2 is the most accurate for ui, as can be seen in 
Figure [TOl For U2 AES-FEM 1 or AES-FEM 2 is the most accurate, and for U 3 
GFD is the most accurate; see Figure ITT] In all these cases, AES-FEM is more 
accurate than FEM. 


6.1.3 Element-Quality Dependence Test 

We test how FEM, AES-FEM, and GFD perform on a series of progressively 
worse meshes. We begin with the most refined mesh from mesh series 2. We 
select 6 of the 523,148 elements and incrementally move one of their nodes 
towards the opposite edge so as to create flatter triangles. We then solve the 
Poisson equation with the polynomial analytic solution ui in (|42l) and record 
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Figure 8 : The errors for 2D convection-diffusion equation on mesh 1 for ui. The 
errors were computed using the Loo norm (left) and the L 2 norm (right). 




Figure 9: The Loo norm errors for the 2D convection-diffusion equation on mesh 
1 for U 2 (left) and U 3 (right). 




Figure 10: The errors for 2D convection-diffusion equation on mesh 2 for ui- 
The errors were computed using the Loo norm (left) and the L 2 norm (right). 
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Figure 11: The Loo norm errors for the 2D convection-diffusion equation on 
mesh 2 for U 2 (left) and U 3 (right). 


the condition numbers of the coefficient and stiffness matrices and the numbers 
of iterations required for the solver to converge. Since the stiffness matrix is the 
same for AES-FEM 1 and AES-FEM 2, the results are just labeled as AES-FEM. 
We use the conjugate gradient method with incomplete Cholesky preconditioner 
for FEM and we use GMRES with incomplete LU preconditioner for AES-FEM 
and GFD. The tolerance for the solvers is 10“® and the drop tolerance for the 
preconditioners is 10“^. As a measure of the mesh quality, we consider the 
cotangent of the minimum angle in the mesh; as the minimum angle tends to 
zero, the cotangent tends towards infinity. Eor very small angles, the cotangent 
of the angle is approximately equal to the reciprocal of the angle. We estimate 
the condition numbers using the MATLAB function condest, which computes 
a lower bound for the 1 -norm condition number. 

The worse the mesh quality, the higher the condition number of the stiffness 
matrix resulting from FEM. In contrast, the condition numbers of the GFD 
coefficient matrix and stiffness matrix from AES-FEM remain almost constant. 
As the condition number for EEM rises, so does the number of iterations re¬ 
quired for the solver to converge, from 102 to 128. The numbers of iterations 
required to solve the equation for AES-EEM and GFD remain constant, at 73 
and 74, respectively. We show the results for 6 meshes. Preconditioned conju¬ 
gate gradient stagnates when trying to solve the FEM linear system from the 
7th mesh, where the minimum angle is approximately 9.1 x 10“® degrees. Solv¬ 
ing the AES-FEM and GED linear systems from the 7th mesh requires the same 
numbers of iterations as the other meshes, 73 and 74 respectively. See Figure 
M for a comparison of the condition numbers and the numbers of iterations. 

The errors for both AES-FEM 1 and AES-FEM 2 rose slightly between the 
3rd and the 4th mesh; the errors were then constant for the rest of the meshes. 
The errors for EEM remained constant and the errors for GFD remained nearly 
constant over the 6 meshes. AES-FEM 1, AES-FEM 2 and GFD converge on 
the 7th mesh in the series with the same errors as on Mesh 6 , whereas for FEM 
the solver stagnates. See Table [T] for specific errors. 
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Table 1: Errors in L 2 norm for FEM, AES-FEM 1, AES-FEM 2, and GFD for 
Ml on a series of meshes with progressively worse mesh element quality. 



FEM 

AES-FEM 1 

AES-FEM 2 

GFD 

Mesh 1 

2.42 X 

2.47 X 10"*^ 

7.82 X 10"'^ 

1.11 X 10-‘> 

Mesh 2 

2.42 X 10“^ 

2.47 X 10"*^ 

7.83 X 10-' 

1.10 X lO-*^ 

Mesh 3 

2.42 X 10“*^ 

2.47 X lO-*^ 

7.82 X 10-^ 

1.10 X 10"^* 

Mesh 4 

2.42 X 10““ 

2.81 X 

1.19 X lO-*^ 

1.10 X 10-« 

Mesh 5 

2.42 X 10-*^ 

2.81 X lO"*^ 

1.19 X lO-^' 

1.10 X 10-^* 

Mesh 6 

2.42 X 10-*^ 

2.81 X lO-*^ 

1.19 X 10-*^ 

1.10 X 10-« 




Figure 12: The condition numbers of the stiffness matrices for FEM and Adap¬ 
tive Extended Stencil (AES)-FEM and the coefficient matrix for generalized fi¬ 
nite difference (GFD) (left) and the numbers of solver iterations (right). Solvers 
used are preconditioned conjugate gradient (PGG) for FEM and preconditioned 
generalized minimal residual (GMRES) for AES-FEM and GFD. 
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Figure 13: Runtimes for a 2D convection-diffusion equation on the most refined 
mesh in series 2. 


6.1.4 Efficiency 

We compare the runtimes for the four methods: AES-FEM 1, AES-FEM 2, 
FEM, and GFD. We consider the convection-diffusion equation on the most 
refined mesh of series 2 with the polynomial analytic solution ui for this runtime 
experiment. We decompose the total time into 4 subcategories: Initialization, 
which includes the time to load the mesh and the time to assign the boundary 
conditions and problem values; Assembly, which includes the time to build the 
stiffness matrix and load vector; Preconditioner, which is the time it takes to 
construction the matrix preconditioner using incomplete LU factorization with 
a drop tolerance of 10“^; and Solver, which is the amount of time for solving 
the preconditioned system using GMRES with a tolerance of 10“®. See Figure 
m for the comparison. The initialization time is minuscule compared to the 
other categories and is not visible in the figure. FEM requires 76 iterations of 
GMRES to converge, AES-FEM 1 and AES-FEM 2 both require 74 iterations, 
and GFD requires 75 iterations. 

One can see that FEM has an advantage when it comes to efficiency on the 
same mesh. In terms of total runtime, it is about 2.1 times faster than AES- 
FEM 1, 2.1 times faster than AES-FEM 2, and 1.8 times faster than GFD. In 
terms of assembly time, FEM is about 6.1 times faster than AES-FEM 1, 7.7 
times faster than AES-FEM 2, and 2.2 times faster than GFD. Comparing the 
assembly time of AES-FEM and GFD, we see that GFD is approximately 3.2 
and 3.5 times faster than AES-FEM 1 and AES-FEM 2, respectively. 

In 2D, assembling the load vector using FEM basis functions (AES-FEM 1) 
saves some time compared to using GLP basis functions (AES-FEM 2). AES- 
FEM 1 offers a savings of approximately 1.1 seconds or, in other words, a 10.3% 
reduction of assembly time and a 3.5% reduction of total time compared to 
AES-FEM 2. We will see in the next section that the efficiency of these two 
methods varies more in 3D. 

However, in terms of error versus runtime, AES-FEM is competitive with, 
and often is more efficient than, the classical FEM with linear basis functions. 
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Figure 14: Loo norm errors versus runtimes for a 2D Poisson equation (left) and 
convection-diffusion equation (right) on mesh series 2. Lower is better. 


In Figure 1141 we compare the Loo norm errors versus runtimes for the four 
methods on mesh series 2 for the Poisson equation and the convection-diffusion 
equation with the analytic solution equal to U 2 . For the Poisson equation, all 
four methods are very similar, with AES-FEM 1 being slightly more efficient for 
finer meshes. For the convection-diffusion equation, AES-FEM 1 and AES-FEM 
2 are approximately the same in terms of efficiency and are the most efficient. 
GFD is also more efficient than FEM. 

6.2 3D Results 

In this section, we present the results from the 3D experiments. We consider 
the Poisson equation and the convection-diffusion equation in 3D. We test three 
problems for each equation on two different series of meshes, each with four 
levels of refinement. The first series of meshes (referred to collectively as “mesh 
series 1”) is created by placing nodes on a regular grid and using MATLAB’s 
Delaunay triangularization to create the elements. The meshes in series 1 range 
from 8x8x8 nodes to 64 x 64 x 64 nodes. The minimum dihedral angle in the 
most refined mesh of series 1 is 35.2 degrees and the maximum dihedral angle 
is 125.2 degrees. The maximum aspect ratio is 4.9, where the aspect ratio of 
a tetrahedron is defined as the ratio of the longest edge length to the smallest 
height. The second series of meshes (referred to collectively as “mesh series 2”) 
is created using TetGen [3]. The number of nodes in mesh series 2 for each 
level of refinement is 509, 4,080, 32,660, and 261,393, which is approximately 
the same as the meshes in series 1. The minimum dihedral angle of the most 
refined mesh in series 2 is 6.7 degrees and the largest dihedral angle is 165.5 
degrees. The largest aspect ratio is 15.2. 
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Figure 15: The errors for 3D Poisson equation on mesh 1 with ui. The errors 
were computed using the Lao norm (left) and the L 2 norm (right). 

6.2.1 Poisson Equation 

We first consider the Poisson equation with Dirichlet boundary conditions on 


the unit cube. That is, 

-V^u = / in D, (47) 

u = g on dLl. (48) 

where fl = [0,1]^. We consider three different analytic solutions listed below. 

ui = 64a: {1 — x) y {1 — y) z {1 — z), (49) 

U 2 = cos( 7 ra;) cos( 7 ry) cos( 7 r 2 :), (50) 


U 3 = -- 5 - - —sinh( 7 rx) cosh( 7 rj/) cosh( 7 rz). (51) 

smh TT cosh tt cosh tt 

The Dirichlet boundary conditions are derived from the analytic solutions. They 
are homogeneous for ui and non-homogeneous for U 2 and U 3 . 

The Loo and L 2 norm errors for the Poisson equation on mesh series 1 for 
ui can be seen in Figure [15] and the Loo norm errors for 1 x 2 and U 3 can be seen 
in Figure [TBl GFD is the most accurate for ui and U 2 - AES-FEM 2 is the most 
accurate for U 3 . 

For mesh series 2, AES-FEM 2 is the most accurate and is close to an order 
of magnitude more accurate than FEM. See Figure |T7| for the Loo and L 2 norm 
errors for ui and Figure |TH| for the Loo norm errors for U 2 and U 3 . 


6.2.2 Convection-Diffusion Equation 

We consider the convection-diffusion equation with Dirichlet boundary condi¬ 
tions on the unit cube, D = [ 0 , 1 ]^. 

—-I- c • Vu = / 
u = g 
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in D, 
on dVL. 


(52) 

(53) 
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Figure 16: The Loo norm errors for the 3D Poisson equation on mesh 1 for U 2 
(left) and M 3 (right). 



Number of Vertices 



Figure 17: The errors for 3D Poisson equation on mesh 2 for ui. The errors 
were computed using the Loo norm (left) and the L 2 norm (right). 




Figure 18: The Loo norm errors for the 3D Poisson equation on mesh 2 for U 2 
(left) and U 3 (right). 
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Figure 19: The errors for 3D convection-diffusion equation on mesh 1 for ui. 
The errors were computed using the Loo norm (left) and the L 2 norm (right). 




Figure 20: The Loo norm errors for the 3D convection-diffusion equation on 
mesh 1 for U 2 (left) and U 3 (right). 

We take c = [1,1,1]^ and we consider the same analytic solutions as in the 
previous section. Again, the Dirichlet boundary conditions are derived from the 
analytic solutions iti, it 2 , and U3. 

The Loo and L 2 norm errors for the 3D convection-diffusion equation on 
mesh series 1 for ui can be seen in Figure fTOl see Figure [50] for the Loo norm 
errors on mesh series 1 for U2 and U3. As with the Poisson equation, GFD is the 
most accurate for ui and U 2 - For U 3 , either GFD or AES-FEM 2 is the most 
accurate depending on the level of refinement. 

On mesh series 2, AES-FEM 2 is the most accurate for ui and U 2 , as can 
be seen in Figure [5T] and the left panel of Figure [551 For U 3 , either GFD or 
AES-FEM 2 is the most accurate, as can be seen in the right panel of Figure [551 
Similar to 2D results, AES-FEM is more accurate than FEM in all these cases. 
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Figure 21: The errors for 3D convection-diffusion equation on mesh 2 for ui. 
The errors were computed using the Too norm (left) and the L 2 norm (right). 




Figure 22: The Loo errors for the 3D convection-diffusion equation on mesh 2 
for U 2 (left) and U 3 (right). 
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6.2.3 Element-Quality Dependence Test 

We test how FEM, AES-FEM, and GFD perforin on a series of meshes with 
progressively worse element shape quality. We begin with the most refined mesh 
from mesh series 1. We select 69 out of the 1,500,282 elements and incrementally 
move one of their nodes towards the opposite side so as to create sliver tetrahe- 
dra. We then solve the Poisson equation with the polynomial analytic solution 
ui in (Si and record the condition numbers of the coefficient and stiffness ma¬ 
trices and the numbers of iterations required for the solver to converge. We 
use the conjugate gradient method with Gauss-Seidel preconditioner for FEM 
and we use GMRES with Gauss-Seidel preconditioner for AES-FEM and GFD. 
We use a tolerance of 10“^ for both solvers. As a measure of the mesh quality, 
we consider the cotangent of the minimum dihedral angle in the mesh; as the 
minimum angle tends to zero, the cotangent tends towards infinity. For very 
small angles, the cotangent of the angle is approximately equal to the reciprocal 
of the angle. We estimate the condition numbers using the MATLAB function 
condest, which computes a lower bound for the 1-norm condition number. 

The worse the mesh quality, the higher the condition number of the stiffness 
matrix resulting from FEM. The condition numbers of the stiffness matrix from 
AES-FEM and the coefficient matrix from GFD remain almost constant. As the 
condition number of the matrix rises so does the number of iterations required 
for the solver to converge. For FEM the number of iterations increases from 69 
for the best mesh to 831 for the most deformed mesh. The numbers of iterations 
for AES-FEM and GFD remain almost constant, increasing from 56 to 59 and 
from 56 to 60, respectively. See Figure [23] for a comparison of the condition 
numbers and iteration counts of the solvers. 

For each of the four methods, the errors were nearly constant over the series 
of meshes. For FEM, the L 2 error on the 1st mesh was 4.36 x 10“'’’ and on the 
6 th mesh, the error was 4.37 x 10“^. For AES-FEM 1, the L 2 error on all the 
meshes was 2.92 x 10“^. For AES-FEM 2, the L 2 error on all the meshes was 
1.30 X 10“^. For GFD, the L 2 error on the 1st mesh was 3.43 x 10“® and on 
the 6th mesh, the error was 3.40 x 10“®. 

6.2.4 Efficiency 

We compare the runtimes of the four methods for solving the convection-diffusion 
equation with the polynomial analytic solution ui on the most refined mesh of 
series 2. As in 2D, the total time is decomposed into 4 subcategories: Ini¬ 
tialization, Assembly, Preconditioner, and Solver. The preconditioner used is 
incomplete LU with a drop tolerance of 10“^. GMRES with a tolerance of 10“® 
is used as the solver. AES-FEM 1 and AES-FEM 2 each require 142 iterations 
to converge, FEM requires 128 iterations, and GFD requires 138 iterations. The 
majority of the time is spent assembling the matrix and solving the system. See 
Figure [Ml for the comparison. 

As in 2D, FEM is the most efficient method on a given mesh. Overall, the 
total runtime of FEM is approximately 1.7 times faster than AES-FEM 1, 1.9 
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Figure 23: The condition numbers of the stiffness matrices for FEM and AES- 
FEM and the coefficient matrix for GFD (left) and the numbers of solver itera¬ 
tions (right). Solvers used are preconditioned CG for FEM and preconditioned 
GMRES for AES-FEM and GFD. 



Figure 24: Runtimes for a 3D convection-diffusion equation on the most refined 
mesh in series 2. 


times faster than AES-FEM 2, and 1.8 times faster than GFD. The assembly of 
FEM is approximately 5.6 times faster than AES-FEM 1, 6.5 times faster than 
AES-FEM 2, and 1.2 times faster than GFD. 

In 3D, the difference of assembling the load vector using FEM basis func¬ 
tions (AES-EEM 1) versus using GLP basis functions (AES-FEM 2) is more 
pronounced than in 2D. The assembly in AES-FEM 1 is 8.3 seconds shorter 
than that of AES-FEM 2. This means the assembly of AES-FEM 1 uses 14.4% 
less time than that of AES-FEM 1 and the total runtime is 8.7% shorter. 

However, similar to 2D, AES-FEM is competitive with, and most of time 
more efficient than, the classical FEM with linear basis functions in terms of 
error versus runtime. Figure [25] shows the Loo norm errors versus runtimes for 
the four methods on mesh series 2 for the Poisson equation and the convection- 
diffusion equation for U 2 - For the Poisson equation, GFD is more efficient on 
coarser meshes and AES-FEM 2 is more efficient for finer meshes. For the 
convection-diffusion equation, GFD is more efficient on smaller meshes and AES- 
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Figure 25: Lao norm errors versus runtimes for a 3D Poisson equation (left) and 
convection-diffusion equation (right) on mesh series 2. Lower is better. 


FEM 2 is more efficient for finer meshes. AES-FEM 1 is also more efficient than 
FEM. 

7 Conclusions and Future Work 

In this paper, we proposed the adaptive extended stencil finite element method, 
which uses generalized Lagrange polynomial basis functions constructed from 
weighted least squares approximations. The method preserves the theoretical 
framework of the classical FEM and the simplicity in imposing essential bound¬ 
ary conditions and integrating the stiffness matrix. We presented the formula¬ 
tion of AES-FEM, showed that the method is consistent, and discussed both 
the local and global stability of the method. We described the implementation, 
including the mesh data structure and the numerical algorithms. We compared 
the accuracy of AES-FEM against the classical EEM with linear basis func¬ 
tions and the quadratic generalized finite difference method for the Poisson and 
convection-diffusion equations in both 2D and 3D. We showed improved accu¬ 
racy and stability of AES-FEM over FEM, and demonstrated that the condition 
number of AES-FEM, and hence the convergence rate of iterative solvers, are 
independent of the element quality of the mesh. Our experiments also showed 
that AES-FEM is more efficient than the classical FEM in terms of error ver¬ 
sus runtime, while having virtually the same sparsity patterns of the stiffness 
matrices. 

As a general method, AES-EEM can use generalized Lagrange polynomial 
basis functions of arbitrary degrees. We only focused on quadratic basis func¬ 
tions in this paper. In future work, we will report higher-order AES-EEM with 
cubic and higher-degree basis functions. The present implementation of AES- 
FEM uses the standard hat functions as the weight functions, which may lead 
to large errors when applied to tangled meshes with inverted elements. We will 
report the resolution of tangled meshes in a future publication. Finally, while 
AES-FEM is efficient in terms of error versus runtime, it is much slower than 
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the classical FEM on a given mesh due to the slower computation of the basis 
functions and the nonsymmetry of the stiffness matrix. The efficiency can be 
improved substantially by leveraging parallelism and efficient multigrid solvers, 
which we will report in the future. 
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